Theta
# ARIMA-Black-Scholes: Semi-Parametric Risk-Neutral Pricing
# T. Moudiki
# 2025-12-04
library(esgtoolkit)
## Loading required package: ggplot2
## Loading required package: gridExtra
## Loading required package: reshape2
## Loading required package: VineCopula
## Loading required package: randtoolbox
## Loading required package: rngWELL
## This is randtoolbox. For an overview, type 'help("randtoolbox")'.
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:zoo':
##
## yearmon, yearqtr
## The following objects are masked from 'package:reshape2':
##
## dcast, melt
##
##
## This is version 1.10.0 of esgtoolkit. Starting with 1.0.0, package renamed as: 'esgtoolkit' (lowercase)
##
##
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(ahead)
# =============================================================================
# 1 - SVJD SIMULATION (Physical Measure)
# =============================================================================
set.seed(123)
n <- 100L
h <- 5
freq <- "daily"
r <- 0.05
maturity <- 5
S0 <- 100
mu <- 0.08
sigma <- 0.04
# Simulate under physical measure with stochastic volatility and jumps
sim_GBM <- esgtoolkit::simdiff(n=n, horizon = h, frequency = freq, x0=S0, theta1 = mu, theta2 = sigma)
sim_SVJD <- esgtoolkit::rsvjd(n=n, r0=mu)
sim_Heston <- esgtoolkit::rsvjd(n=n, r0=mu,
lambda = 0,
mu_J = 0,
sigma_J = 0)
cat("Simulation dimensions:\n")
## Simulation dimensions:
cat("Start:", start(sim_SVJD), "\n")
## Start: 0 1
cat("End:", end(sim_SVJD), "\n")
## End: 5 1
cat("Paths:", ncol(sim_SVJD), "\n")
## Paths: 100
cat("Time steps:", nrow(sim_SVJD), "\n\n")
## Time steps: 1261
par(mfrow=c(1, 3))
# Plot historical (physical measure) paths
esgtoolkit::esgplotbands(sim_GBM, main="GBM Paths under the Physical Measure",
xlab="Time",
ylab="Stock prices")
esgtoolkit::esgplotbands(sim_SVJD, main="SVJD Paths under the Physical Measure",
xlab="Time",
ylab="Stock prices")
esgtoolkit::esgplotbands(sim_Heston, main="Heston Paths under the Physical Measure",
xlab="Time",
ylab="Stock prices")
# Summary statistics
cat("Physical measure statistics (GBM):\n")
## Physical measure statistics (GBM):
terminal_prices_physical_GBM <- sim_GBM[nrow(sim_GBM), ]
cat("Mean terminal price:", mean(terminal_prices_physical_GBM), "\n")
## Mean terminal price: 150.3905
cat("Std terminal price:", sd(terminal_prices_physical_GBM), "\n")
## Std terminal price: 13.40266
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
cat("Physical measure statistics (SVJD):\n")
## Physical measure statistics (SVJD):
terminal_prices_physical_SVJD <- sim_SVJD[nrow(sim_SVJD), ]
cat("Mean terminal price:", mean(terminal_prices_physical_SVJD), "\n")
## Mean terminal price: 150.5358
cat("Std terminal price:", sd(terminal_prices_physical_SVJD), "\n")
## Std terminal price: 18.01574
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
cat("Physical measure statistics (Heston):\n")
## Physical measure statistics (Heston):
terminal_prices_physical_Heston <- sim_Heston[nrow(sim_Heston), ]
cat("Mean terminal price:", mean(terminal_prices_physical_Heston), "\n")
## Mean terminal price: 151.4444
cat("Std terminal price:", sd(terminal_prices_physical_Heston), "\n")
## Std terminal price: 16.65026
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
# =============================================================================
# 2 - COMPUTE DISCOUNTED PRICES (Transform to Martingale Domain)
# =============================================================================
discounted_prices_GBM <- esgtoolkit::esgdiscountfactor(r=r, X=sim_GBM)
discounted_prices_SVJD <- esgtoolkit::esgdiscountfactor(r=r, X=sim_SVJD)
discounted_prices_Heston <- esgtoolkit::esgdiscountfactor(r=r, X=sim_Heston)
martingale_diff_GBM <- discounted_prices_GBM - S0
martingale_diff_SVJD <- discounted_prices_SVJD - S0
martingale_diff_Heston <- discounted_prices_Heston - S0
cat("Martingale differences dimensions (GBM):", dim(martingale_diff_GBM), "\n")
## Martingale differences dimensions (GBM): 1261 100
cat("Mean martingale diff (should be ≠0 under P):\n")
## Mean martingale diff (should be ≠0 under P):
print(t.test(rowMeans(martingale_diff_GBM)))
##
## One Sample t-test
##
## data: rowMeans(martingale_diff_GBM)
## t = 60.9, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 8.356654 8.912984
## sample estimates:
## mean of x
## 8.634819
cat("\nMartingale differences dimensions (SVJD):", dim(martingale_diff_SVJD), "\n")
##
## Martingale differences dimensions (SVJD): 1261 100
cat("Mean martingale diff (should be ≠0 under P):\n")
## Mean martingale diff (should be ≠0 under P):
print(t.test(rowMeans(martingale_diff_SVJD)))
##
## One Sample t-test
##
## data: rowMeans(martingale_diff_SVJD)
## t = 57.891, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 8.102235 8.670648
## sample estimates:
## mean of x
## 8.386441
cat("\nMartingale differences dimensions (Heston):", dim(martingale_diff_Heston), "\n")
##
## Martingale differences dimensions (Heston): 1261 100
cat("Mean martingale diff (should be ≠0 under P):\n")
## Mean martingale diff (should be ≠0 under P):
print(t.test(rowMeans(martingale_diff_Heston)))
##
## One Sample t-test
##
## data: rowMeans(martingale_diff_Heston)
## t = 58.357, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 8.384831 8.968209
## sample estimates:
## mean of x
## 8.67652
# =============================================================================
# 3 - VISUALIZE RISK PREMIUM
# =============================================================================
par(mfrow=c(2,2))
matplot(discounted_prices_GBM, type='l', col=rgb(0,0,1,0.3),
main="Discounted Stock Prices (Physical Measure - GBM)",
ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)
matplot(discounted_prices_SVJD, type='l', col=rgb(0,0,1,0.3),
main="Discounted Stock Prices (Physical Measure - SVJD)",
ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)
matplot(discounted_prices_Heston, type='l', col=rgb(0,0,1,0.3),
main="Discounted Stock Prices (Physical Measure - Heston)",
ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
mean_disc_path_GBM <- rowMeans(discounted_prices_GBM)
times_plot <- as.numeric(time(discounted_prices_GBM))
plot(times_plot, mean_disc_path_GBM, type='l', lwd=2, col='blue',
main="Risk Premium in Discounted Prices (GBM)",
xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
col=c('blue','red','green'), lty=c(1,2,3), lwd=2)
mean_disc_path_SVJD <- rowMeans(discounted_prices_SVJD)
times_plot <- as.numeric(time(discounted_prices_SVJD))
plot(times_plot, mean_disc_path_SVJD, type='l', lwd=2, col='blue',
main="Risk Premium in Discounted Prices (SVJD)",
xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
col=c('blue','red','green'), lty=c(1,2,3), lwd=2)
mean_disc_path_Heston <- rowMeans(discounted_prices_Heston)
times_plot <- as.numeric(time(discounted_prices_Heston))
plot(times_plot, mean_disc_path_Heston, type='l', lwd=2, col='blue',
main="Risk Premium in Discounted Prices (Heston)",
xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
col=c('blue','red','green'), lty=c(1,2,3), lwd=2)
# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================
n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)
martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)
# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()
arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()
arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()
# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_GBM[, i])
fit <- forecast::thetaf(y)
arima_models_GBM[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_GBM[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_SVJD[, i])
fit <- forecast::thetaf(y)
arima_models_SVJD[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_SVJD[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_Heston[, i])
fit <- forecast::thetaf(y)
arima_models_Heston[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_Heston[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}
# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM),
function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
##
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.5035183
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.96
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD),
function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
##
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.4609734
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.84
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston),
function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
##
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
main="Box-Ljung P-values (GBM)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_SVJD, breaks=20, col='lightblue',
main="Box-Ljung P-values (SVJD)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_Heston, breaks=20, col='lightcoral',
main="Box-Ljung P-values (Heston)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================
bootstrap_independent <- function(x, R = 1000, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
n <- length(x)
# Each column is an independent bootstrap sample
replicate(R, sample(x, size = n, replace = TRUE))
}
cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
##
##
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20 # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)
# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()
# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i],
R = n_sim_per_path)
fit <- arima_models_GBM[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_GBM[[i]] <- S_tilde_price
}
# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i],
R = n_sim_per_path)
fit <- arima_models_SVJD[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_SVJD[[i]] <- S_tilde_price
}
# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i],
R = n_sim_per_path)
fit <- arima_models_Heston[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_Heston[[i]] <- S_tilde_price
}
# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM),
frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD),
frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston),
frequency = frequency(sim_Heston))
# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM,
main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD,
main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston,
main="Risk-Neutral Paths - Heston")
# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (GBM)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (SVJD)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (Heston)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================
cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
##
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)
cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 128.3881
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.01439695
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_GBM - capitalized_stock_price
## t = -0.051539, df = 1999, p-value = 0.9589
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.5622233 0.5334294
## sample estimates:
## mean of x
## -0.01439695
cat("\nSVJD Risk-Neutral Verification:\n")
##
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 128.1574
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.2451161
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_SVJD - capitalized_stock_price
## t = -0.61866, df = 1999, p-value = 0.5362
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.0221290 0.5318968
## sample estimates:
## mean of x
## -0.2451161
cat("\nHeston Risk-Neutral Verification:\n")
##
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 128.5378
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: 0.1352577
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_Heston - capitalized_stock_price
## t = 0.40195, df = 1999, p-value = 0.6878
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.5246852 0.7952007
## sample estimates:
## mean of x
## 0.1352577
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================
cat("\n=== OPTION PRICING ===\n\n")
##
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {
d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
d2 <- d1 - sigma * sqrt(T)
call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
put <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)
list(call = call, put = put)
}
strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)
# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
n_strikes <- length(strikes)
call_prices <- numeric(n_strikes)
bs_call_prices <- numeric(n_strikes)
put_prices <- numeric(n_strikes)
bs_put_prices <- numeric(n_strikes)
call_se <- numeric(n_strikes)
put_se <- numeric(n_strikes)
for (k in 1:n_strikes) {
K <- strikes[k]
call_payoffs <- pmax(terminal_prices - K, 0)
call_prices[k] <- mean(call_payoffs) * discount_factor
call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
put_payoffs <- pmax(K - terminal_prices, 0)
put_prices[k] <- mean(put_payoffs) * discount_factor
put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
bs_call_prices[k] <- bs_prices$call
bs_put_prices[k] <- bs_prices$put
}
list(call_prices = call_prices, put_prices = put_prices,
bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
call_se = call_se, put_se = put_se)
}
# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
| 37.6847250 |
0.0000000 |
37.6959374 |
0.0000001 |
0.2175495 |
0.0000000 |
| 29.8980644 |
0.0013472 |
29.9079897 |
0.0000602 |
0.2174527 |
0.0013472 |
| 22.1452991 |
0.0365898 |
22.1260239 |
0.0061022 |
0.2154791 |
0.0092665 |
| 14.6435444 |
0.3228429 |
14.4726589 |
0.1407450 |
0.2036358 |
0.0336206 |
| 8.0122810 |
1.4795873 |
7.6644048 |
1.1204988 |
0.1710185 |
0.0788649 |
| 3.2699519 |
4.5252661 |
3.0014135 |
4.2455153 |
0.1165797 |
0.1375942 |
| 0.9066440 |
9.9499659 |
0.8280730 |
9.8601826 |
0.0606846 |
0.1860634 |
| 0.1685717 |
16.9999015 |
0.1608817 |
16.9809992 |
0.0227524 |
0.2096258 |
| 0.0121927 |
24.6315303 |
0.0226357 |
24.6307609 |
0.0053122 |
0.2167928 |
kableExtra::kable(as.data.frame(options_Heston))
| 37.8028720 |
0.0015958 |
37.6959374 |
0.0000001 |
0.2619526 |
0.0015958 |
| 30.0399788 |
0.0267105 |
29.9079897 |
0.0000602 |
0.2603824 |
0.0089609 |
| 22.3576823 |
0.1324218 |
22.1260239 |
0.0061022 |
0.2552281 |
0.0240528 |
| 15.0348491 |
0.5975964 |
14.4726589 |
0.1407450 |
0.2383996 |
0.0534645 |
| 8.6676523 |
2.0184075 |
7.6644048 |
1.1204988 |
0.2021327 |
0.1015914 |
| 4.0766962 |
5.2154591 |
3.0014135 |
4.2455153 |
0.1478793 |
0.1598164 |
| 1.5391986 |
10.4659694 |
0.8280730 |
9.8601826 |
0.0911697 |
0.2103637 |
| 0.4459361 |
17.1607147 |
0.1608817 |
16.9809992 |
0.0486331 |
0.2421993 |
| 0.1148506 |
24.6176371 |
0.0226357 |
24.6307609 |
0.0241640 |
0.2554788 |
kableExtra::kable(as.data.frame(options_SVJD))
| 37.5880705 |
0.0830298 |
37.6959374 |
0.0000001 |
0.3025649 |
0.0233010 |
| 29.9091959 |
0.1921630 |
29.9079897 |
0.0000602 |
0.2964688 |
0.0395838 |
| 22.3743079 |
0.4452828 |
22.1260239 |
0.0061022 |
0.2853585 |
0.0617545 |
| 15.2673169 |
1.1262996 |
14.4726589 |
0.1407450 |
0.2627206 |
0.0947876 |
| 9.1376458 |
2.7846363 |
7.6644048 |
1.1204988 |
0.2229804 |
0.1415378 |
| 4.5976251 |
6.0326235 |
3.0014135 |
4.2455153 |
0.1689933 |
0.1972375 |
| 1.9332399 |
11.1562461 |
0.8280730 |
9.8601826 |
0.1140386 |
0.2462269 |
| 0.7097549 |
17.7207690 |
0.1608817 |
16.9809992 |
0.0705917 |
0.2786469 |
| 0.2313729 |
25.0303948 |
0.0226357 |
24.6307609 |
0.0434949 |
0.2958461 |
ETS
# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================
n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)
martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)
# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()
arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()
arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()
# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_GBM[, i])
fit <- forecast::ets(y)
arima_models_GBM[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_GBM[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_SVJD[, i])
fit <- forecast::ets(y)
arima_models_SVJD[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_SVJD[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_Heston[, i])
fit <- forecast::thetaf(y)
arima_models_Heston[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_Heston[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}
# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM),
function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
##
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.5048272
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.96
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD),
function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
##
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.4611386
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.84
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston),
function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
##
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
main="Box-Ljung P-values (GBM)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_SVJD, breaks=20, col='lightblue',
main="Box-Ljung P-values (SVJD)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_Heston, breaks=20, col='lightcoral',
main="Box-Ljung P-values (Heston)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================
bootstrap_independent <- function(x, R = 1000, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
n <- length(x)
# Each column is an independent bootstrap sample
replicate(R, sample(x, size = n, replace = TRUE))
}
cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
##
##
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20 # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)
# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()
# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i],
R = n_sim_per_path)
fit <- arima_models_GBM[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_GBM[[i]] <- S_tilde_price
}
# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i],
R = n_sim_per_path)
fit <- arima_models_SVJD[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_SVJD[[i]] <- S_tilde_price
}
# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i],
R = n_sim_per_path)
fit <- arima_models_Heston[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_Heston[[i]] <- S_tilde_price
}
# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM),
frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD),
frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston),
frequency = frequency(sim_Heston))
# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM,
main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD,
main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston,
main="Risk-Neutral Paths - Heston")
# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (GBM)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (SVJD)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (Heston)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================
cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
##
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)
cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 127.924
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.4785094
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_GBM - capitalized_stock_price
## t = -1.7301, df = 1999, p-value = 0.08377
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.0209334 0.0639145
## sample estimates:
## mean of x
## -0.4785094
cat("\nSVJD Risk-Neutral Verification:\n")
##
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 127.8367
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.5658013
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_SVJD - capitalized_stock_price
## t = -1.3747, df = 1999, p-value = 0.1694
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.3729490 0.2413465
## sample estimates:
## mean of x
## -0.5658013
cat("\nHeston Risk-Neutral Verification:\n")
##
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 128.3782
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: -0.02433147
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_Heston - capitalized_stock_price
## t = -0.072939, df = 1999, p-value = 0.9419
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.6785473 0.6298844
## sample estimates:
## mean of x
## -0.02433147
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================
cat("\n=== OPTION PRICING ===\n\n")
##
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {
d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
d2 <- d1 - sigma * sqrt(T)
call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
put <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)
list(call = call, put = put)
}
strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)
# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
n_strikes <- length(strikes)
call_prices <- numeric(n_strikes)
bs_call_prices <- numeric(n_strikes)
put_prices <- numeric(n_strikes)
bs_put_prices <- numeric(n_strikes)
call_se <- numeric(n_strikes)
put_se <- numeric(n_strikes)
for (k in 1:n_strikes) {
K <- strikes[k]
call_payoffs <- pmax(terminal_prices - K, 0)
call_prices[k] <- mean(call_payoffs) * discount_factor
call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
put_payoffs <- pmax(K - terminal_prices, 0)
put_prices[k] <- mean(put_payoffs) * discount_factor
put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
bs_call_prices[k] <- bs_prices$call
bs_put_prices[k] <- bs_prices$put
}
list(call_prices = call_prices, put_prices = put_prices,
bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
call_se = call_se, put_se = put_se)
}
# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
| 37.3233395 |
0.0000657 |
37.6959374 |
0.0000001 |
0.2153985 |
0.0000657 |
| 29.5427541 |
0.0074881 |
29.9079897 |
0.0000602 |
0.2148244 |
0.0053024 |
| 21.7860973 |
0.0388391 |
22.1260239 |
0.0061022 |
0.2130421 |
0.0128624 |
| 14.2477281 |
0.2884778 |
14.4726589 |
0.1407450 |
0.2029315 |
0.0332499 |
| 7.6572043 |
1.4859618 |
7.6644048 |
1.1204988 |
0.1702183 |
0.0777218 |
| 3.0857449 |
4.7025103 |
3.0014135 |
4.2455153 |
0.1154527 |
0.1362044 |
| 0.8569653 |
10.2617385 |
0.8280730 |
9.8601826 |
0.0608434 |
0.1841160 |
| 0.1686439 |
17.3614249 |
0.1608817 |
16.9809992 |
0.0251250 |
0.2069742 |
| 0.0200318 |
25.0008207 |
0.0226357 |
24.6307609 |
0.0091880 |
0.2140408 |
kableExtra::kable(as.data.frame(options_Heston))
| 37.6822414 |
0.0052534 |
37.6959374 |
0.0000001 |
0.2593967 |
0.0032356 |
| 29.9177530 |
0.0287729 |
29.9079897 |
0.0000602 |
0.2578997 |
0.0110308 |
| 22.2663566 |
0.1653843 |
22.1260239 |
0.0061022 |
0.2510846 |
0.0277011 |
| 14.9579035 |
0.6449390 |
14.4726589 |
0.1407450 |
0.2334187 |
0.0579578 |
| 8.5691693 |
2.0442127 |
7.6644048 |
1.1204988 |
0.1971864 |
0.1052936 |
| 3.9707040 |
5.2337551 |
3.0014135 |
4.2455153 |
0.1421274 |
0.1627970 |
| 1.4344219 |
10.4854809 |
0.8280730 |
9.8601826 |
0.0850850 |
0.2126207 |
| 0.3989463 |
17.2380131 |
0.1608817 |
16.9809992 |
0.0416221 |
0.2426567 |
| 0.0714887 |
24.6985633 |
0.0226357 |
24.6307609 |
0.0170685 |
0.2558071 |
kableExtra::kable(as.data.frame(options_SVJD))
| 37.3936091 |
0.1383182 |
37.6959374 |
0.0000001 |
0.3098695 |
0.0393111 |
| 29.7336475 |
0.2663644 |
29.9079897 |
0.0000602 |
0.3029494 |
0.0551093 |
| 22.2189875 |
0.5397123 |
22.1260239 |
0.0061022 |
0.2912961 |
0.0767336 |
| 15.1512978 |
1.2600304 |
14.4726589 |
0.1407450 |
0.2679371 |
0.1088503 |
| 9.0474757 |
2.9442161 |
7.6644048 |
1.1204988 |
0.2286590 |
0.1542833 |
| 4.5472247 |
6.2319730 |
3.0014135 |
4.2455153 |
0.1762522 |
0.2081396 |
| 1.9296726 |
11.4024287 |
0.8280730 |
9.8601826 |
0.1245120 |
0.2553862 |
| 0.7570340 |
18.0177979 |
0.1608817 |
16.9809992 |
0.0854171 |
0.2860005 |
| 0.3383181 |
25.3870899 |
0.0226357 |
24.6307609 |
0.0595079 |
0.3010065 |
ARIMA
# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================
n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)
martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)
# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()
arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()
arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()
# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_GBM[, i])
fit <- forecast::auto.arima(y)
arima_models_GBM[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_GBM[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_SVJD[, i])
fit <- forecast::auto.arima(y)
arima_models_SVJD[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_SVJD[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}
# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
y <- as.numeric(martingale_increments_Heston[, i])
fit <- forecast::thetaf(y)
arima_models_Heston[[i]] <- fit
res <- as.numeric(residuals(fit))
arima_residuals_Heston[, i] <- res
centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}
# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM),
function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
##
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.6636805
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.95
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD),
function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
##
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.6953426
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.98
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston),
function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
##
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
main="Box-Ljung P-values (GBM)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_SVJD, breaks=20, col='lightblue',
main="Box-Ljung P-values (SVJD)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
hist(pvalues_Heston, breaks=20, col='lightcoral',
main="Box-Ljung P-values (Heston)",
xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================
bootstrap_independent <- function(x, R = 1000, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
n <- length(x)
# Each column is an independent bootstrap sample
replicate(R, sample(x, size = n, replace = TRUE))
}
cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
##
##
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20 # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)
# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()
# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i],
R = n_sim_per_path)
fit <- arima_models_GBM[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_GBM[[i]] <- S_tilde_price
}
# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i],
R = n_sim_per_path)
fit <- arima_models_SVJD[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_SVJD[[i]] <- S_tilde_price
}
# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i],
R = n_sim_per_path)
fit <- arima_models_Heston[[i]]
fitted_increments <- as.numeric(fitted(fit))
discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
discounted_path[1, ] <- S0
increments <- matrix(scale(fitted_increments, center = TRUE,
scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) +
resampled_residuals[1:(n_periods - 1), ]
discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
S_tilde_price <- discounted_path * discount_factor
all_S_tilde_Heston[[i]] <- S_tilde_price
}
# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM),
frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD),
frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston),
frequency = frequency(sim_Heston))
# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM,
main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD,
main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston,
main="Risk-Neutral Paths - Heston")
# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (GBM)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (SVJD)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)],
type='l', col=rgb(0,0,1,0.1),
main="Risk-Neutral Paths (Heston)",
xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================
cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
##
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)
cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 127.8758
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.5267259
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_GBM - capitalized_stock_price
## t = -1.825, df = 1999, p-value = 0.06816
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.09275998 0.03930816
## sample estimates:
## mean of x
## -0.5267259
cat("\nSVJD Risk-Neutral Verification:\n")
##
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 128.2571
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.1454399
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_SVJD - capitalized_stock_price
## t = -0.35471, df = 1999, p-value = 0.7228
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.9495529 0.6586730
## sample estimates:
## mean of x
## -0.1454399
cat("\nHeston Risk-Neutral Verification:\n")
##
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 127.9524
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: -0.4501411
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
##
## One Sample t-test
##
## data: terminal_prices_rn_Heston - capitalized_stock_price
## t = -1.3619, df = 1999, p-value = 0.1734
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1.0983259 0.1980436
## sample estimates:
## mean of x
## -0.4501411
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (GBM)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (SVJD)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
main="Terminal Prices: Physical (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)
hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
main="Terminal Prices: Risk-Neutral (Heston)",
xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================
cat("\n=== OPTION PRICING ===\n\n")
##
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {
d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
d2 <- d1 - sigma * sqrt(T)
call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
put <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)
list(call = call, put = put)
}
strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)
# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
n_strikes <- length(strikes)
call_prices <- numeric(n_strikes)
bs_call_prices <- numeric(n_strikes)
put_prices <- numeric(n_strikes)
bs_put_prices <- numeric(n_strikes)
call_se <- numeric(n_strikes)
put_se <- numeric(n_strikes)
for (k in 1:n_strikes) {
K <- strikes[k]
call_payoffs <- pmax(terminal_prices - K, 0)
call_prices[k] <- mean(call_payoffs) * discount_factor
call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
put_payoffs <- pmax(K - terminal_prices, 0)
put_prices[k] <- mean(put_payoffs) * discount_factor
put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
bs_call_prices[k] <- bs_prices$call
bs_put_prices[k] <- bs_prices$put
}
list(call_prices = call_prices, put_prices = put_prices,
bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
call_se = call_se, put_se = put_se)
}
# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
| 37.2870565 |
0.0013337 |
37.6959374 |
0.0000001 |
0.2246654 |
0.0013337 |
| 29.5063372 |
0.0086222 |
29.9079897 |
0.0000602 |
0.2241400 |
0.0057295 |
| 21.7615753 |
0.0518681 |
22.1260239 |
0.0061022 |
0.2217762 |
0.0145639 |
| 14.2625991 |
0.3408998 |
14.4726589 |
0.1407450 |
0.2104212 |
0.0372080 |
| 7.7729434 |
1.6392520 |
7.6644048 |
1.1204988 |
0.1759057 |
0.0826745 |
| 3.2759164 |
4.9302327 |
3.0014135 |
4.2455153 |
0.1193116 |
0.1418863 |
| 0.9644973 |
10.4068215 |
0.8280730 |
9.8601826 |
0.0623878 |
0.1912890 |
| 0.1787918 |
17.4091238 |
0.1608817 |
16.9809992 |
0.0242475 |
0.2163885 |
| 0.0157002 |
25.0340400 |
0.0226357 |
24.6307609 |
0.0059463 |
0.2238247 |
kableExtra::kable(as.data.frame(options_Heston))
| 37.3495392 |
0.0041721 |
37.6959374 |
0.0000001 |
0.2570826 |
0.0030129 |
| 29.5933710 |
0.0360118 |
29.9079897 |
0.0000602 |
0.2550569 |
0.0116692 |
| 21.9374713 |
0.1681198 |
22.1260239 |
0.0061022 |
0.2485708 |
0.0279110 |
| 14.6507009 |
0.6693573 |
14.4726589 |
0.1407450 |
0.2302646 |
0.0585080 |
| 8.3516001 |
2.1582644 |
7.6644048 |
1.1204988 |
0.1917575 |
0.1070118 |
| 3.7789938 |
5.3736659 |
3.0014135 |
4.2455153 |
0.1363246 |
0.1653930 |
| 1.2901727 |
10.6728526 |
0.8280730 |
9.8601826 |
0.0804251 |
0.2145030 |
| 0.3303330 |
17.5010208 |
0.1608817 |
16.9809992 |
0.0398569 |
0.2426599 |
| 0.0673177 |
25.0260132 |
0.0226357 |
24.6307609 |
0.0186772 |
0.2534207 |
kableExtra::kable(as.data.frame(options_SVJD))
| 37.7713083 |
0.1886397 |
37.6959374 |
0.0000001 |
0.3023061 |
0.0587409 |
| 30.1244079 |
0.3297471 |
29.9079897 |
0.0000602 |
0.2943756 |
0.0732991 |
| 22.6284991 |
0.6218461 |
22.1260239 |
0.0061022 |
0.2812424 |
0.0937678 |
| 15.4824587 |
1.2638136 |
14.4726589 |
0.1407450 |
0.2592663 |
0.1231765 |
| 9.2623385 |
2.8317012 |
7.6644048 |
1.1204988 |
0.2205829 |
0.1645298 |
| 4.6409015 |
5.9982720 |
3.0014135 |
4.2455153 |
0.1668160 |
0.2151494 |
| 1.8977393 |
11.0431176 |
0.8280730 |
9.8601826 |
0.1118326 |
0.2617140 |
| 0.6901214 |
17.6235076 |
0.1608817 |
16.9809992 |
0.0688966 |
0.2916385 |
| 0.2333527 |
24.9547468 |
0.0226357 |
24.6307609 |
0.0402776 |
0.3074406 |